<!DOCTYPE html>
<html lang="zh-cn">
<head>
    <title>线性微分方程组</title>
    <meta charset="utf-8" />
    <link rel="stylesheet" type="text/css" href="../../css/note.css" />
</head>
<body>

<p>
	本章研究方程
	<span class="formula">
		`x' = A(t) x + f(t)`,
	</span>
	`A` 为矩阵, `x, x', f` 为向量.
</p>

<h2>常数变易公式</h2>

<p>	若已知齐方程 `x' = A(t)x` 的基解矩阵为 `Phi(t)`, 于是
	`Phi'(t) = A(t) Phi(t)`. 可设
	`varphi(t) = Phi(t) c(t)` 是非齐方程的解, 代入得
	<span class="formula">
		`Phi(t) c'(t) = f(t)`.
	</span>
	因为 `Phi(t)` 可逆, 所以令
	<span class="formula">
		`c(t) = int_(t_0)^t Phi^-1(s) f(s) "d"s`,
	</span>
	就得到一个满足初值条件 `varphi(t_0) = 0` 的解. 一般地, 满足初值条件
	`varphi(t_0) = eta` 的解为
	<span class="formula">
		`varphi(t) = Phi(t) (Phi^-1(t_0) eta
					+ int_(t_0)^t Phi^-1(s) f(s) "d"s)`.
	</span>
</p>

<h2>n 阶非齐次线性方程的常数变易公式</h2>

<p>	因为
	<span class="formula">
		`sum_(i=0)^n a_i(t) x^((n-i)) = f(t)`, `quad a_0(t) = 1`,
		`x^((i))(t_0) = eta_(i+1)`, `i = 0, 1, cdots, n-1`
		<span class="label" id="for-nth-linear-nonhomogen">(5-1)</span>
	</span>
	等价于
	<span class="formula">
		`x' = Ax + vec(f)(t)`, `quad x(t_0) = eta`,
	</span>
	其中 `x = (x_1, x_2, cdots, x_n)^T = (x, x', cdots, x^((n-1)))^T`,
	`eta = (eta_1, eta_2, cdots, eta_n)^T`,
	`A = [
			0, 1, 0, cdots, 0;
			0, 0, 1, cdots, 0;
			vdots, vdots, vdots, , vdots;
			0, 0, 0, cdots, 1;
			-a_n(t), -a_(n-1)(t), -a_(n-2)(t), cdots, -a_1(t);
	]`,
	`vec(f)(t) = [0; 0; vdots; 0; f(t)]`. <br/>
	则由线性方程组的常数变易公式知, 若 <a class="ref"
		href="#for-nth-linear-nonhomogen"></a> 对应的齐方程的基本解组为
	`x_1, x_2, cdots, x_n`, 则 <a class="ref"
		href="#for-nth-linear-nonhomogen"></a> 有特解
	<span class="formula">
		`varphi(t) = sum_(k=1)^n x_k(t)
					 int_(t_0)^t (W_k(s))/(W(s)) f(s) "d"s`,
	</span>
	其中
	<span class="formula">
		`W(s) = |
			x_1(s), x_2(s), cdots, x_n(s);
			x_1'(s), x_2'(s), cdots, x_n'(s);
			vdots, vdots, , vdots;
			x_1^((n-1))(s), x_2^((n-1))(s), cdots, x_n^((n-1))(s);
		|`,<br/>
		`W_k(s)` 是 `W(s)` 中 `x_k^((n-1))(s)` 的代数余子式.
	</span>
	特别当 `n = 2`, 上式化为
	<span class="formula">
		`varphi(t) = int_(t_0)^t
		|x_1(s), x_2(s); x_1(t), x_2(t)| /
		|x_1(s), x_2(s); x_1'(s), x_2'(s)| f(s) ds`.
	</span>
</p>

<h2>常系数线性方程组</h2>

<p>
	<span class="formula">
		`x' = Ax`
	</span>
</p>

<p>	此方程有基解矩阵 `exp At`. 显然 `exp A0 = E`.</p>

<p>	若 `lambda` 是 `A` 的特征值, `bm v` 为对应的特征向量, 则
	`"e"^(lambda t) bm v` 是一个解. 特别地, 若 `A` 有 `n`
	个线性无关的特征向量, 对应的特征值为 `lambda_1, cdots, lambda_n`,
	则方程有基解矩阵 `["e"^(lambda_1 t) bm v_1, cdots, "e"^(lambda_n t)
	bm v_n]`.
</p>

<p>	若已知一个基解矩阵 `Phi(t)`, 则存在一常数矩阵 `C`, 使得
	`exp At = Phi(t) C`, 再由 `exp A0 = E` 知 `exp At = Phi(t) Phi^-1(0)`.
</p>

<h3>利用根子空间求矩阵指数</h3>

<p>
  矩阵指数
  <span class="formula">
    `exp bm A = sum_(k ge 0) bm A^k/(k!)`
  </span>
  是无穷级数, 直接计算比较困难.  但如果存在某个 `k_0`, 使得 `bm A = bm O`,
  `AA k ge k_0`, 则只需计算有限项.  根子空间就具有这样的性质. 具体来说,
  设 `lambda` 是 `bm A` 的一个特征值, 重数为 `r_lambda`,
  <span class="formula">
    `V_lambda = "Ker"(bm A-lambda bm E)^(r_lambda)`
    `:= {bm v: (bm A-lambda bm E)^(r_lambda) bm v = bb 0}`
  </span>
  是对应的根子空间 (特别 `r_lambda = 1` 时, 根子空间就是特征子空间),
  则对任意 `bm v in V_lambda` 有
  <span class="formula">
    `(bm A-lambda bm E)^k bm v = bb 0`, `AA k ge r_lambda`.
  </span>
  根据线性代数知识, `n` 维线性空间在根子空间上具有直和分解, 即任意
  `n` 维向量 `bm x` 可以唯一写成
  <span class="formula">
    `bm x = sum_lambda bm x_lambda`, `quad bm x_lambda in V_lambda`.
  </span>
  上式表示对 `bm A` 的不同特征值求和. 从而
  <span class="formula">
    `bm x exp bm A`
    `= sum_lambda bm x_lambda exp bm A`
    `= sum_lambda "e"^lambda bm x_lambda exp(bm A-lambda bm E)`
    `= sum_lambda "e"^lambda bm x_lambda sum_(k lt r_lambda) (bm A-lambda
    bm E)^k // k!`.
  </span>
  特别取 `bm x = bm epsi_j`, 就得到矩阵 `exp bm A` 的第 `j` 列.
</p>

<p class="example">
  求 `exp bm A t`, `bm A = [cos a, -sin a; sin a, cos a]`.
  推论: `exp[, -a; a,] = [cos a, -sin a; sin a, cos a]`.<br/>
</p>

<p class="solution">
  `bm A` 的特征值是 `lambda_(1,2) = cos a +- "i" sin a`,
  先设 `lambda_1 = lambda_2 = cos a`, 则 `bm A = lambda bm E`,
  <span class="formula">
    `exp(bm A - lambda bm E) = exp bm O = bm E`.
  </span>
  从而
  <span class="formula">
    `[1;0] exp bm A t = [1;0] "e"^(lambda t)`,<br/>
    `[0;1] exp bm A t = [0;1] "e"^(lambda t)`,<br/>
    `exp bm A t = "e"^(t cos a) bm E`.
  </span>
  再设 `lambda_1 != lambda_2`, 对应的特征向量为
  <span class="formula">
    `bm v_1 = ["i";1]`, `bm v_2 = [1;"i"]`.
  </span>
  将单位向量分解到特征子空间,
  <span class="formula">
    `[1;0] = 1/2 [1;"i"] - "i"/2 ["i";1]`,<br/>
    `[0;1] = -"i"/2 [1;"i"] + 1/2 ["i";1]`.
  </span>
  由于特征值都是一重的, 指数的无穷级数中实际上只有一项:
  <span class="formula">
    `exp(bm A - lambda bm E) = bm E`.
  </span>
  从而
  <span class="formula">
    `[1;0] exp bm A t`
    `= 1/2 [1;"i"] "e"^(lambda_2 t) - "i"/2 ["i";1] "e"^(lambda_1 t)`
    `= 1/2 ["e"^(lambda_2 t) + "e"^(lambda_1 t); "i"("e"^(lambda_2 t)
    - "e"^(lambda_1 t))]`
    `= "e"^(t cos a) [cos(t sin a); sin(t sin a)]`.
  </span>
  同理可求得第二列, 于是
  <span class="formula">
    `exp bm A t`
    `= "e"^(t cos a) [cos(t sin a), -sin(t sin a);
    sin(t sin a), cos(t sin a)]`.
  </span>
  易知上式对重根的情形也成立.
</p>

<p class="example">
	用线性微分方程组的知识解 Lanchester 方程 (见第 4 章例题).
</p>

<p class="solution">
	矩阵 `A = [0, -beta^2; -alpha^2, 0]` 的特征值为 `+-alphabeta`,
	对应的特征向量为 `(beta, &mp; alpha)^T`.
	代入公式计算,
	<span class="formula">
		`exp At (1;0) = 1/(2beta)("e"^(alpha beta t) (beta;-alpha)
		+ "e"^(-alpha beta t) (beta;alpha))`,<br/>
		`exp At (0;1) = 1/(2alpha)("e"^(alpha beta t) (-beta;alpha)
		+ "e"^(-alpha beta t) (beta;alpha))`.
	</span>
	所以基解矩阵
	<span class="formula">
		`Phi(t) = exp At = [
			cosh(alpha beta t), -beta/alpha sinh(alpha beta t);
			-alpha/beta sinh(alpha beta t), cosh(alpha beta t);
	]`.</span>
	而
	<span class="formula">
		`Phi^(-1)(0) = [
			cosh(alpha beta t), beta/alpha sinh(alpha beta t);
			alpha/beta sinh(alpha beta t), cosh(alpha beta t);
		](0) = E`.
	</span>
	故 Lanchester 方程满足初值条件 `(A, B)^T(0) = (A_0, B_0)^T` 的解为
	<span class="formula">
		`Phi(t)Phi^(-1)(0) (A_0, B_0)^T = Phi(t)(A_0, B_0)^T`.
	</span>
</p>

<p class="example">
  寻找级数 `sum_(n ge 0) x^(3n)/((3n)!)` 的闭形式.
</p>

<p class="solution">
  记 `x_i = sum_(n ge 0) t^(3n+i)/((3n+i)!)`, `i = 0, 1, 2`, 于是
  `x_0' = x_2`, `x_2' = x_1`, `x_1' = x_0`. 写成向量
  `bm x = (x_0, x_1, x_2)^T`, 有
  <span class="formula">
    `bm x' = bm (A x)`, `quad bm A = [0, 0, 1; 1, 0, 0; 0, 1, 0]`.
  </span>
  容易看出 `exp bm A t = bm E x_0(t) + bm A x_1(t) + bm A^2 x_2(t)`,
  但这只会使我们回到原点, 因此我们寻找它的另一种表示.
  `bm A` 的特征多项式为 `lambda^3 - 1`, 三个根为 `1, omega, omega^2`,
  对应的特征向量为
  <span class="formula">
    `bm v_1 = [1; 1; 1]`,
    `bm v_2 = [1; omega^2; omega]`,
    `bm v_3 = [1; omega; omega^2]`.
  </span>
  注意到 `x_0(0) = 1, x_1(0) = x_2(0) = 0`, 因此方程组的初值条件恰为
  单位向量 `bm eta = (1, 0, 0)^T`. 将这个向量分解为
  <span class="formula">
    `bm eta = 1/3 (bm v_1 + bm v_2 + bm v_3)`,
  </span>
  从而解得
  <span class="formula">
    `bm x`
    `= bm eta exp bm A t`
    `= 1/3 (bm v_1 "e"^t + bm v_2 "e"^(omega t) + bm v_3 "e"^(omega^2 t))`
    `= 1/3 [
      "e"^t + 2 "e"^(t/2) cos((sqrt 3)/2 t);
      "e"^t + 2 "e"^(t/2) cos((sqrt 3)/2 t - (2pi)/3);
      "e"^t + 2 "e"^(t/2) cos((sqrt 3)/2 t + (2pi)/3)
    ]`.
  </span>
</p>

<script src="../../js/note.js?type=math"></script>
</body>
</html>
